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THERMOPLASTIC  ANALYSIS  USING  THE  FINITE  ELEMENT 


CODE  AGGIE  I 

David  H.  Allen*  and  Walter  E.  Haisler** 

ABSTRACT 

The  authors  have  previously  proposed  a  theoretical  model  for 
predicting  structural  response  of  elastic-plastic-creep  materials 
subjected  to  variable  temperature  loadings.  This  model  util  Ires 
the  combined  Isotropic-kinematic  hardening  rule  In  conjunction 
with  the  classical  rate  Independent  plasticity  theory.  The 
resulting  constitutive  law  has  been  Implemented  Into  the  finite 
strain  finite  element  code  AGGIE  I.  The  authors  have  found  it 
necessary  to  alter  the  formulation  previously  presented  In  order 
to  obtain  a  concise  numerical  formulation.  That  modification  is 
presented  In  this  paper.  In  addition,  a  section  Is  included 
which  outlines  supplementary  information  necessary  for  computer 
program  Input,  followed  by  a  computational  procedure  which  Is 
Intended  to  aid  programmers  who  wish  to  use  this  model.  Finally, 
the  results  of  several  sample  cases  are  presented. 

INTRODUCTION 

In  our  publication  entitled  "The  Application  of  Thermal  and  Creep 
Effects  to  the  Combined  Isotropic-Kinematic  Hardening  Model  for  Inelastic 
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Structural  Analysis  by  the  Finite  Element  Method"1,  we  presented  an 
infinitesimal  strain  constitutive  law  intended  for  use  with  elastic- 
plastic  materials  subjected  to  thermal  loadings.  The  model  accounted 
for  time  dependent  behavior  through  an  uncoupled  creep  term.  The  result¬ 
ing  law  was  utilized  in  a  finite  element  formulation  which  uses  a  total 
Lagrangian  formulation,  and  the  resulting  equations  of  motion  were 
presented  at  that  time. 

The  impetus  for  this  new  law  resulted  from  a  survey  of  the  literature 

on  temperature  dependent  constitutive  laws.  It  was  determined  that  most 

finite  element  programs  available  which  account  for  thermal  effects  utilize 

either  isotropic  or  kinematic  hardening  and  thus  may  be  inadequate  for 

modelling  the  Bauschinger  effect  during  cyclic  loading.  This  statement 

5 

has  been  previously  verified  for  isothermal  cases  by  Haisler  et  al . 

We  found  subsequent  to  our  first  publication  that  a  similar  model 
had  been  previously  proposed  by  Yamada1^’1^.  However,  his  model  contains 
significant  variations  from  that  proposed  by  us.  Specifically,  he  does 
not  account  for  elastic  modulus  changes  due  to  temperature,  and  his 
model  is  constructed  within  an  Eulerian  framework.  More  importantly, 
his  determination  of  the  hardening  modulus  departs  significantly  from 
that  presented  herein.  This  latter  point  will  be  elaborated  upon  later 
in  this  report. 

During  implementation  of  our  model  into  the  finite  element  code 
AGGIE  I6,  we  found  it  necessary  to  make  a  significant  modification 
to  the  model  in  order  to  obtain  computational  results.  This  alteration 
as  well  as  a  short  review  of  the  model  will  be  presented  in  the  following 
section. 

An  added  purpose  of  this  publication  is  to  present  the  supplementary 
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Information  necessary  to  input  the  model  to  a  computer  program.  Thus, 
we  have  Included  a  section  detailing  these  additions,  as  well  as  a 
computational  procedure  Intended  to  clarify  the  model  for  those 
who  wish  to  pursue  the  topic  In  greater  depth. 

In  the  final  section  the  results  of  several  sample  problems  are 
presented  as  a  verification  of  the  completed  computer  program. 

REVIEW  AND  MODIFICATION  OF  THE  THEORY 

Recall  that  in  the  incremental  theory  of  plasticity,  the  work¬ 
hardening  rule  Is  defined  by  a  yield  function  in  stress  space,  which 
for  temperature  dependent  combined  isotropic-kinematic  hardening  is 

F(Sij  '  c,1j)  3  k?  (  ;d‘P’T)* 

where 

=  second  Piola-Kirchhoff  stress  tensor, 

«  coordinates  of  the  yield  surface  center  in  stress  space, 

_p 

de^  »  uniaxial  plastic  strain  increment, 

T  *  temperature, 

n 

and  / dt  represents  the  plastic  strain  history  dependence  in  the  yield 
function.  Note  that  the  explicit  temperature  dependence  on  the  right 
hand  side  of  equation  (1)  precludes  kinematic  hardening  due  to  temp¬ 
erature  changes.  A  schematic  representation  of  equation  (1)  is  shown 
in  Figure  1. 

Differentiating  equation  (1)  gives  the  following  consistency 
condition  during  plastic  loading: 
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(2) 


where  the 
which  can 


ts,j  dSu  •  d“i.t  - 2k  ^ djP-  2klr 111  •  °> 


3F  i  -  a,.) 

term  a?  represents  — r*^ - evaluated  at 

”  i j  351j 

be  seen  to  be  equivalent  to,|£  x  . 

t3bij  "  °i j  ' 


Since  during  neutral  loading  the  plastic  strain  increment  and  da^j  are 
zero,  it  is  apparent  that  a  statement  governing  loading  is 


Ht/Su  -  2klr  "T  *>• 

In  addition  to  the  statement  of  consistency  during  loading  given  by 
(2),  the  following  associated  flow  rule  is  employed: 


(3) 


*  dX 


3F 

SS-  • 
U 


The  stress  at  any  time  is  assumed  to  be  given  by 


(4) 


S1j  '  Dijmn  ^n  '  Emn  "  Emn  Emn^ 


where  is  the  elastic  constitutive  tensor  at  time  t. 
of  equation  (5)  was  shown  in  our  previous  paper^  to  give 


(5) 

Incrementation 


dS1J  ’  DV)m  <«»,  ■  dC  -  dEL-  dO 
+  ^mn  '  ^mn  ^mn^* 

where 


(6) 


^ijrnn  *  e^*s^c  modulus  tensor  at  time  t+At, 

dD1Jmn  *  change  in  the  elastic  constitutive  tensor  (due  to  a 
J  temperature  change)  during  time  step  At, 
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E^  ■  total  strain  tensor, 

p 

E^  *  plastic  strain  tensor, 
r 

E^  *  creep  strain  tensor, 

E^  •  thermal  strain  tensor, 

and  all  superscripts  t  Indicate  measurement  at  time  t.  The  superscript 
t+At  on  the  elastic  constitutive  tensor  Is  necessary  due  to  the  fact  that 
the  time  step  is  finite  rather  than  Infinitesimal. 

The  final  expression  required  to  complete  the  constitutive  law 
Is  the  hardening  rule.  According  to  Ziegler's  modification  a  tensorially 
correct  statement  is 

da1j  "  dv,<S1j  "  a1j)*  (7) 

where  u  Is  a  scalar  to  be  determined  by  the  consistency  condition 
[equation  (2)]. 

The  purpose  of  the  above  equations  is  to  obtain  a  mapping  giving 
the  total  stress  Increment  In  terms  of  the  total  strain  Increment  and 
known  quantities.  To  accomplish  this  it  is  necessary  to  solve  equations 
(2),  (4) ,  (6)  and  (7)  to  obtain  the  stress  Increment  tensor.  Unfortunately, 

p 

these  comprise  nineteen  equations  in  twenty  unknowns:  dS^.dE^,  do^.dx, 

and  dv>.  One  would  normally  consider  this  problem  underspecified,  but 
It  will  turn  out  that  dx  is  an  artificial  parameter  which  need  not  be 
determined  In  order  to  obtain  a  solution.  Incidentally,  It  Is  possible 
to  obtain  dx  after  the  constitutive  law  Is  solved  by  determining  the  plastic 
strain  increment. 

The  constitutive  law  Is  thus  obtained  by  substituting  equation 
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(4)  Into  equations  (2)  and  (6)  and  solving  for  the  stress  Increment 
tensor.  The  resulting  equation  Is 


nt+At  3F  3F  nt+At 
Ijvw  turnn 


tu  +  2  k—  4^ 
x  ,-P  <fiT 

at 


3F 

3^ 


tu 


+  2k 


3k  dT 

si  ar 


+  D 


t+At  3F 
pqrs  3lpq 


x 


("wC  <> 


^JULdD*“- 


- IB - 

OL  3k  de  .  3k  dT  .  nt+At  3F  3F 

!?  ^  ^ 31  W"  Spq  ^rs. 


+  dD 


(8) 


The  above  constitutive  law  differs  from  that  obtained  by  Yamada16,17 
In  that  the  term  containing  dT  Is  written  as  a  separate  portion  of  the 
stress  Increment  In  the  latter's  derivation.  It  will  be  shown,  however, 
that  the  two  are  mathematically  equivalent.  It  should  also  be  noted 
that  the  above  derivation  differs  slightly  from  that  proposed  In  our 
previous  paper.  We  have  chosen  this  new  form  due  to  Its  clarity, 
although  either  formulation  Is  mathematically  correct.  Unfortunately, 

In  the  form  presented  In  equation  (8)  It  Is  not  possible  to  determine 
the  stress  Increment  due  to  the  occurrence  of  several  undetermined 
parameters  on  the  right  hand  side.  Yamada  has  accounted  for  these 
terms  by  Introducing  the  plastic  work  rate.  We  have  chosen  a  slightly 
different  method  for  determining  the  effect  of  these  unknowns. 
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As  presented  In  a  previous  paper  by  us  ,  we  assume  that  there  exists 
a  scaler  parameter  c,  called  the  hardening  modulus,  which  when  multiplied 
by  the  plastic  strain  Increment  and  subtracted  from  the  stress  Increment 
will  be  parallel  to  a  tangent  to  the  yield  surface.  Mathematically,  this 
may  be  stated  as 


(dS^  -  cdEP)  |£  -  0  . 

ij  U  (g) 

It  can  be  seen  from  an  examination  of  the  consistency  condition  [equation 
(2)],  that  in  order  for  the  above  statement  to  be  valid,  the  following 
must  hold: 


cdEij  =  daij  +  2k^fp  dl?  +  2kll  dT 


If  one  employs  the  normality  condition  in  equation  (10)  it  is  found 


3F  3F  _  3F  datu  .  ?.3k  deP  _.3k  dT 

Xn  "  *♦..  “35T  +  2k7>  3T  +  2k3l  dx 

pq  pq  tu  3e  , 

which  when  substituted  into  the  constitutive  relation  [equation  (8)] 


results  in 


0t+At  |£  rjt+At 

uijmn  3\w  3btu  Utumn 
-3F  3F  .  nt+At  3F  3T” 
Ws  3*  3?  c 

pq  pq  ™  pq  rs 


<dE™,  -  <„  -  <„> 


t+At  3F  3F  \ 

D1jvw  ^vw  *tu  d^tumn  \  ,Ft  _  FPt  FCt  _  FTt, 
.aF IF  nt+At  dF  aF  I  ^Lmn  "  Lmn  "  Lmn  "  Swr 
Ca5nfl  a?  +  Dpqrs  3?  a?  / 

pq  pq  ^  pq  rs  / 


♦  dW^  -  C  -  &  -  C>  . 
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The  simplicity  of  the  above  formulation  can  be  seen  when  one  uses 
the  normality  condition  In  conjunction  with  equation  (9)  to  obtain 


c 


dS 


3F 

*^1j 

ir~dE 


u 

r 

u 


03) 


Thus,  for  Isothermal  loading  It  can  be  seen  that  the  hardening  modulus 
is  simply  two  thirds  the  instantaneous  slope  of  the  uniaxial  stress- 
plastic  strain  diagram.  However,  for  non- Isothermal  loadings  this  Is  not 
the  case.  To  see  this  note  that  since  the  uniaxial  stress  Is  a  function 
of  both  the  plastic  strain  history  and  temperature,  an  Increment  of 
uniaxial  stress  is  given  by 


A  —  ^  j“P  ,  30  Jr 

do  *  — p  de  +  dl 
dc 


04) 


Combining  equations  (13)  and  (14)  gives 


c 


2 

7 


(H*  + 


do  dT  \ 
dc  , 


(15) 


where  H'  is  the  Instantaneous  slope  of  the  stress-plastic  strain  diagram. 
Substituting  the  above  equation  Into  the  constitutive  relation  [equation 
(12)]  ,  then  rearranging  and  employing  the  normality  condition,  gives 


dStj  ■  C?WdEmn  *  d4  -  <>  *  dP 


U 


(16) 


where 


it+At 


Dt+At  I?  If  nt+4t 

uljvw  _  vw  3btu  utumn 

'ijmn  Ijmn  "  2  u,  3F 3F  .  nt+At  3F  sF 

3  P^S  ^ 

pq  pq  r'l  pq  « 
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dP 
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dT 


(18) 


and  the  terms  H’  and  ’y  niay  be  determined  from  uniaxial  stress-strain 
data  at  time  t  as  described  later.  Thus,  it  can  be  seen  that  c  is 
superscripted  at  time  t  since  the  hardening  parameters  are  detennined 
there.  Due  to  its  length  the  derivation  of  the  above  constitutive 
equation  is  not  covered  in  detail  here,  but  may  be  obtained  from  the 
authors  on  request. 

The  translation  of  the  yield  surface  in  stress  space  may  now  be 
obtained  by  substituting  equation  (7)  into  equation  (2)  and  solving 
for  dp.  The  resulting  relation  is 


O  L-  3  K  JT  _i  —  P 


du  = 


If  dS  -  2k^  dT  -  2kp,  £ 
JJ _  3t 


(S 


mn 


mn 


3F 

3$ 


mn  .  (19) 

Equation  (19)  thus  guarantees  that  the  state  of  stress  will  remain 
consistent  with  the  yield  surface  during  loading  even  under  non- isothermal 

conditions. 

To  summarize,  then,  the  constitutive  law  is  obtained  by  solving. 
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in  the  following  order,  equations  (18),  (17),  (16),  (19)  and  (7). 

An  application  of  the  constitutive  law  to  the  principle  of  virtual 
work  within  a  total  Lagrangian  framework  yields  the  following  matrix 
algebra  equations: 

M  I:ri+  ([Kl]t  +  [KNL]t)  {AUc}-  M  .  (20 

where  the  quantities  in  the  above  equation  are  as  defined  in  reference  1. 
Solving  the  above  equations  and  employing  equilibrium  iteration  gives 
the  displacement  increment  for  a  given  load  step. 
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APPLICATION  OF  THE  THEORY  TO  THE 
COMPUTER  CODE  AGGIE  I 

Although  the  above  model  Is  mathematically  correct,  certain 
additions  must  be  made  in  order  to  obtain  a  computational  form.  In 
this  section  are  presented  those  additional  requirements  which  were 
necessary  to  input  the  model  to  the  code  AGGIE  I.  It  should  be 
emphasized,  however,  that  these  additions  are  not  essential  to  the  model 
itself,  and,  as  such,  variations  can  be  made  without  affecting  the 
validity  of  the  model. 

A.  The  Yield  Function 

A  function  describing  the  onset  of  permanent  deformation  must  be 
utilized  to  describe  equation  (1).  The  yield  function  utilized  in 
AGGIE  I  is  the  Von  Mises  yield  criterion  for  Isotropic  materials, 
which  is  given  by 

1  2  ,2  2  2 
F(tJi )  =  j  [(oi  -  02)  +  (02  "  °3 )  +  (03  -  <*i)  ]  =  k  ,  (21) 

where  olto 2  and  a3  are  principal  stresses  and  k  is  the  current  uniaxial 

yield  surface  size. 

B.  Creep  Strain  Calculation 

In  this  theory  the  creep  strain  is  considered  to  be  an  initial  strain 
which  can  be  uncoupled  from  the  Instantaneous  plastic  strain.  There  are 
currently  three  methods  of  calculating  the  creep  strain  in  the  program 
AGGIE  I.  One  Is  the  use  of  empirical  equations  from  microphenomenological 
materials  science  literature.  The  second  consists  of  interpolating  a 
set  of  creep  strain  versus  time  curves  at  varying  stress  levels.  The 
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final  method  makes  use  of  viscoelasticity  theory.  Because  of  the 
complexity  of  these  calculations,  further  discussion  Is  omitted  from 
this  paper.  The  reader  is  referred  to  references  (13)  and  (14)  for 
a  more  detailed  discussion. 

C.  Determination  of  Thermal  Strain 

For  the  form  of  the  constitutive  law  given  by  equation  (15)  it  can 
be  shown  that  the  thermal  strain  for  isotropic  materials  can  be 
approximated  by 


T 

ij 


where  ay  is  a  temperature  dependent  material  properly,  TR  Is  the 


(22) 


reference  temperature  and  6^  is  the  Kroneker  delta.  In  metals  literature 
it  is  customary  to  define 


so  that 


=  TT^T 


(23) 


Eij  =  6ija(T  '  V  .  (24) 

It  can  be  seen  that  whenever  ay  is  independent  of  temperature  ay  and  a 


are  identical.  Since  both  ay  and  a  are  reported  in  various  materials 
handbooks,  the  option  to  use  either  is  specified  in  the  code  AGGIE  I. 
When  ay  is  employed  the  thermal  strain  is  calculated  by  integrating 


equation  (22),  and  when  a  is  chosen  the  thermal  strain  is  determined 
by  solving  equation  (24). 
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D.  Material  Input  Data 

In  order  to  facilitate  easy  Input,  the  materials  data  are  Input  as 
piecewise  linear  curves.  The  necessary  data  are  the  uniaxial  stress- 
strain  curves  for  various  temperatures,  a  coefficient  of  thermal  expansion 
versus  temperature  curve,  and  a  Poisson's  ratio  versus  temperature  curve. 
These  data  are  shown  schematically  In  Figure  2. 

E.  Calculation  of  Equivalent  Uniaxial  Yield  Surface  Size  and  Gradients. 

We  have  found  It  convenient  to  construct  a  set  of  k  vs.  TP  curves 


as  shown  In  Figure  3  from  the  Input  stress-strain  data  shown  In  Figure  2. 
This  Is  accomplished  by  a  transformation  on  both  the  uniaxial  stress  and 
strain.  The  i**1  abscissa  may  be  obtained  from 


-P 

E1  "Scl 


°xi  /Ej 


(25) 


where  E^  Is  the  uniaxial  elastic  modulus  of  a  stress-strain  curve  at  1^. 

1  L 

The  1U  ordinate  is  determined  using 


k1  “  °xj  +  0  *  (ox1  “°xj*  ,  (26) 

where  is  the  Initial  yield  stress  at  temperature  T^  and  p  Is  the 

ratio  of  Isotropic  to  kinematic  hardening,  p  may  be  obtained  from  a 
uniaxial  cyclic  loading  test.  Util izlng  equations  (25)  and  (26)  will 
allow  one  to  construct  the  data  shown  in  Figure  3. 

In  order  to  obtain  the  current  value  of  k  It  Is  first  necessary 
to  determine  the  uniaxial  plastic  strain,  which  Is  given  by 


Thus,  for  the  current  temperature  and  uniaxial  plastic  strain  level  It 
Is  possible  to  obtain  k,  as  shown  In  Figure  4.  A  linear  Interpolation 
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gives 


<VT> 


k  '  kH  ■  (kH  ■  kl) 


(28) 


where  T  Is  the  current  temperature.  In  addition.  It  Is  possible  to 
obtain  the  yield  surface  temperature  gradient,  which  Is 

ak  _  kH  '  *L 


(29) 


_ 

It  should  be  noted  that  may  be  obtained  In  a  similar  fashion  from  a 

_D 

stress  vs  uniaxial  plastic  strain  curve.  Also,  one  may  estimate  3k/ 3e  by 


3k 

3e 


(§V[(^)h'0?)JtTO 


(30) 


and  H'  may  be  determined  similarly  from  a  stress  vs  uniaxial  plastic 
strain  curve. 

F.  Elastic  Strain  Increment 

The  elastic  strain  increment  is  obtained  by  solving  equation  (6): 
<n  •  0I?;«  [dS1j  -  dD1J«,  -  C  -  C  -  .  (31 


where  OI^j^  represents  the  elastic  compliance  tensor  at  time  t+At. 


COMPUTATIONAL  PROCEDURE 

Using  the  above  modifications  and  extensions  of  the  theoretical 
model  It  Is  possible  to  establish  a  computational  procedure  which  when 
programmed  will  completely  define  the  constitutive  law. 

Below  Is  the  computational  procedure  for  non- Isothermal  combined 
Isotropic-kinematic  hardening  used  by  AGGIE  I.  The  procedure  Is 
used  In  an  equilibrium  Iteration  loop  and  Is  therefore  performed  several 


-14- 


times  for  each  load  step.  The  first  time  through  the  procedure  during 
each  load  step  the  stresses  are  calculated  for  the  previous  load  step 
In  A,  through  0.  The  materials  data  are  then  updated  to  the  current 
step  In  P.  through  X.  These  data  are  then  used  In  another  routine 
to  establish  the  total  strain  increment.  On  each  succeeding  Iteration 
the  stresses  are  calculated  for  the  current  step  In  A.  through  0.  and 
steps  P.  through  X.  are  not  performed.  The  program  variables  used  In 
this  outline  are  defined  In  a  table  following  the  procedure. 

A.  If  this  Is  the  first  load  step  and  the  first  Iteration  on  that  load 
step  (KSTEP.  EQ.  1.  AND.  1C0UNT.  NE.  3)  calculate  the  initial 
temperature  change. 

B.  Obtain  change  In  elastic  constitutive  tensor  due  to  temperature 
change  during  current  step. 

1.  Interpolate  Input  materials  data  to  obtain  elastic  modulus, 
Poisson's  ratio,  and  coefficient  of  thermal  expansion  for  temperature 
at  start  of  current  step. 

2.  Construct  elastic  constitutive  tensor  for  temperature  at  start 
of  current  step. 

3.  Interpolate  Input  materials  data  to  obtain  elastic  modulus 
and  Poisson's  ratio  for  temperature  at  end  of  current  step. 

4.  Construct  elastic  constitutive  tensor  for  temperature  at  end 
of  current  step. 

5.  Subtract  elastic  constitutive  tensor  at  start  of  step  from 
elastic  constitutive  tensor  at  end  of  step  to  obtain  change  in 
elastic  constitutive  tensor  during  current  step. 

C.  Calculate  the  total  strain  Increment  for  current  step. 

D.  Calculate  thermal  strain  Increment  for  current  step. 
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1.  If  IALPHA  ■  1  the  coefficient  of  thermal  expansion  Is  a 
function  of  temperature  [equation  (22)]. 

2.  If  IALPHA  •  0  the  coefficient  of  thermal  expansion  Is  not  a 
function  of  temperature  [equation  (24)]. 

E.  Obtain  Initial  estimate  of  elastic  strain  Increment. 

1.  Subtract  creep  strain  Increment  from  total  strain  Increment. 

2.  Subtract  thermal  strain  Increment  from  total  strain  Increment. 

3.  Initialize  elastic  strain  Increment  to  elastic  plus  plastic 
strain  Increment. 

F.  Calculate  a  trial  stress  Increment  assuming  the  total  strain  Increment 
to  be  purely  elastic  [equation  (6)]. 

G.  Calculate  trial  stress  state. 

H.  Check  whether  trial  stress  state  causes  yielding  [equation  (1)]. 

If  F>k2  and  IPEL  ■  1  go  to  I. 

If  F>k2  and  IPEL  -  2  go  to  K. 

If  F  *  k2  and  IPEL  ■  1  set  IPEL  ■  2  and  go  to  N. 

If  F  «  k2  and  IPEL  -  2  go  to  K. 

2 

If  F<k  go  to  N. 

I.  The  material  was  behaving  elastically  at  the  beginning  of  this  time 
step,  but  will  yield  during  this  strain  Increment.  Determine  the 
part  of  the  strain  Increment  taken  elastically.  Set  IPEL  ■  2  and 
IELAS  -  2. 

J.  Update  stress  state  to  yield  surface. 

K.  The  material  Is  now  behaving  plastically. 

1.  Determine  the  number  of  subincrements,  M. 

2.  Determine  the  strain  subincrement. 
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3.  Determine  the  temperature  subincrement. 

4.  Interpolate  materials  data  to  obtain  elastic  modulus  and 
Poisson's  ratio  at  end  of  subincrement. 

5.  Construct  elastic  constitutive  tensor  at  end  of  subincrement. 

6.  Subtract  elastic  constitutive  tensor  at  start  of  subincrement 
from  elastic  constitutive  tensor  at  end  of  subincrement  to  obtain 
change  in  elastic  constitutive  tensor  during  subincrement. 

L.  Calculate  the  stress  tensor  minus  the  yield  surface  center. 

M.  Obtain  the  plastic  stress  increment  by  looping  through  this  section 
M  times. 

1.  Update  elastic  constitutive  tensor  to  end  of  subincrement. 

2.  Calculate  a  F/aS^j. 

-P  -P 

3.  Calculate  H',  3o/aT,  ak/ac  and  3k/ aT  for  current  value  of  e  . 

4.  Calculate  effective  modulus  tensor  C.jmn  [equation  (17)]. 

5.  Calculate  dP.j  for  current  subincrement  [equation  (18)]. 

6.  Calculate  stress  increment  for  current  subincrement  [equation  (16)]. 

7.  Calculate  the  elastic  strain  subincrement  [equation  (31)]. 

8.  Subtract  the  elastic  strain  subincrement  from  the  total  strain 
subincrement  to  obtain  the  plastic  strain  subincrement. 

9.  Update  total  elastic  strain  increment  for  current  step. 

10.  Update  equivalent  plastic  strain  [equation  (27)]. 

11.  Determine  the  size  of  the  subsequent  yield  surface  [equation  (28)]. 

12.  Calculate  the  new  position  of  the  yield  surface  center  [equations 
(19)  and  (7)]. 

13.  Update  temperature  to  end  of  sub  increment. 

14.  Calculate  stress  tensor  minus  yield  surface  center. 

N.  Update  stresses,  strains  and  elastic  strains  to  end  of  step. 
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O.  Check  for  equilibrium  Iteration.  If  equilibrium  Iterating,  the 
procedure  is  complete.  If  not,  continue. 

P.  Calculate  the  temperature  change  and  temperature  increment  for 
this  step. 

Q.  Obtain  the  change  in  elastic  constitutive  tensor  due  to  temperature 
change  during  current  step. 

1.  Interpolate  input  materials  data  to  obtain  elastic  modulus,  Poisson's 
ratio,  and  coefficient  of  thermal  expansion  for  temperature  at  start 

of  current  step. 

2.  Construct  elastic  constitutive  tensor  for  temperature  at  start 
of  current  step. 

3.  Interpolate  input  materials  data  to  obtain  elastic  modulus  and 
Poisson's  ratio  for  temperature  at  end  of  current  step. 

4.  Construct  elastic  constitutive  tensor  for  temperature  at  end 
of  current  step. 

5.  Subtract  elastic  constitutive  tensor  at  start  of  step  from 
elastic  constitutive  tensor  at  end  of  step  to  obtain  change  in  elastic 
constitutive  tensor  during  current  step. 

6.  Update  elastic  constitutive  tensor  to  end  of  current  step. 

R.  Calculate  thermal  strain  increment  for  current  step. 

1.  I f I ALPHA  *1  the  coefficient  of  thermal  expansion  is  a  function 
of  temperature  [equation  (22)]. 

2.  If  IAl.PHA  *  0  the  coefficient  of  thermal  expansion  is  not  a 
function  of  temperature  [equation  (24)]. 

S.  Calculate  creep  strain  increment  for  current  step. 

T.  Check  for  yielding: 

If  IPEL  »  1  go  to  V. 
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If  IPEL  •  2  go  to  U. 

U.  Calculate  effective  constitutive  tensor  for  plastic  action. 

1.  Find  a  F/JS^. 

2.  Determine  the  stress  Increment  for  the  current  step. 

3.  Find  H',  3o/3 T,  dk/3r**,  and  3k/3T. 

4.  Determine  effective  constitutive  tensor  [equation  (17)]. 

V.  Calculate  new  yield  surface  site.  If  IPFL  ■  2  go  to  X. 

W.  Set  effective  constitutive  tensor  equal  to  elastic  constitutive  tensor 

for  elastic  action. 

X.  Find  dP^  for  current  step  [equation  (18)]. 

PROGRAM  VARIABLES 
KSTEP  -  load  step  number 
ICOUNT  -  equilibrium  iteration  flag 
EQ.3  -  equilibrium  iteration 
NE.3  -  not  equilibrium  Iteration 
IFLAS  -  first  yield  control  flag 
EQ.l  -  material  has  not  yielded 

EQ.2  -  material  has  yielded 

IALPHA  -  thermal  strain  type  flag 
EQ.l  -  constant  alpha 

EQ.2  -  variable  alpha 

IPEL  -  elastic-plastic  control  flag 
EQ.l  -  material  Is  elastic  on  previous  step 

EQ.2  -  material  is  plastic  on  previous  step 
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APPLICATION  OF  THE  MODEL  TO  SPECIAL  CASES 


In  order  to  provide  a  better  understanding  of  the  model  several 
sample  cases  are  considered  In  this  section.  The  purpose  of  these 
cases  Is  both  to  verify  the  accuracy  of  the  constitutive  law 
and  to  provide  several  economical  test  cases  which  may  be  used 
to  debug  the  model  when  It  Is  Input  to  a  finite  element  code. 

The  calculations  are  simplified  by  modelling  the  response  of  uniaxial 
bars  under  small  strain  conditions  and  subjected  to  quasi-static 
loadings.  The  stress  tensor  Is  therefore  replaced  with  the 

uniaxial  stress  a,  and  the  strain  tensor  E^  Is  replaced  with  the 


uniaxial  strain  c. 


Also,  the  elastic  constitutive  tensor  D^mn 


Is 


replaced  with  the  elastic  modulus  E,  and  the  yield  surface  translation 


tensor 


becomes  «. 


For  simplicity,  there  is  assumed  to  be  no 


creep. 

The  stress  Increment  is  first  predicted  by  assuming  zero  plastic 
strain  and  solving  equation  (6): 


do  -  Et+At(de-dtT)  +  (Et+At-  E*)  (ct-«Pt-rTt)  .  (3 2) 

If  the  loading  condition  is  satisfied  [equation  (1)  or  (3)].  then  the 
stress  Increment  Is  updated  by  using  equations  (16)  through  (18)  which 
reduce  to 
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rt+At  jF  jf 
t  —  at 

Jo  JO 

u<  ilE  ?f.  +  rt+At'  jF  jF 

Jo  Jo  Jo  Jo, 


1/  t  Pt 

U  -t  - 


.Tt) 


JF  r t+At  JF  Jo  \ 

3o  J" _ Jp_  JT 1 

>  2 t  +  rt+At  JF  JF  I 

Jo  Jo  Jo  Jo  / 


+  dE(et-cPt-eTt)  +1 


dT  . 


(33) 


The  yield  surface  translation  scalar  [equation  (19)]  reduces  to 


dp  = 


2(o-a)do-2k|£  dT  -  2k~p  deP 

_ _ _ _  jcr _ 

(o-aT2  (o-u) 


(34) 


and  the  yield  surface  translation  increment  tensor  is  simplified  to 

d«  =  du(o-a)  (35) 

By  replacing  oj  in  equation  (21)  with  o,  the  uniaxial  stress,  it  can  be 

seen  that  the  yield  surface  gradient  is  ^  =  2  (o  -  «).  Finally,  the 

elastic  strain  increment  is  determined  from  equation  (31)  to  be 


dc 


E  _  do-dE(ct-cPt-eTt) 

rt+At 


(36) 


A.  Case  Number  1  -  Elastic  Unloading 

The  first  case  involves  an  axial  bar  which  is  loaded  isothermal ly 
to  some  state  of  stress  outside  the  initial  yield  surface  and  then 
slowly  cooled,  as  shown  in  Figure  5.  For  simplicity,  hardening  is 
assumed  to  be  isotropic. 

It  is  appropriate  at  this  point  to  recall  that  due  to  the  nonlinear 
nature  of  the  problem,  the  solution  is  obtained  through  an  iterative 
procedure.  As  described  in  references  7  and  15,  that  method  assumes 
a  strain  increment  and  iterates  to  solution.  It  has  been  shown  to 
converge  for  small  load  increments.  Since  it  is  not  the  intention 
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of  this  paper  to  demonstrate  the  equilibrium  iteration  procedure 
the  method  used  here  will  be  to  assume  a  reasonable  strain  increment 
based  on  engineering  judgement  and  use  that  prediction  to  obtain 
the  stress  and  strain  increments  for  the  load  step.  For  the  conditions 
shown  in  Figure  5,  a  reasonable  assumption  for  the  strain  increment  is 


SL  ♦  £ 


Et  j-t+At 


+  ,iAT 


where  AT  is  understood  to  be  a  negative  quantity.  From  equation  (32) 
the  stress  increment  is  predicted  to  be 


da  =  E 


►At  /  0*  a4 _ \ 

\  E1  Et+At/ 


+  (Et+At-  Et)(3ty  -  1.5  cy)  =  0. 


Since  the  predicted  stress  increment  is  zero  and  the  yield  surface 
size  will  increase  due  to  the  temperature  decrease,  unloading  is 
predicted.  Thus,  there  is  no  plastic  straining  during  the  load  step 
and  the  stress  increment  is  precisely  that  given  by  equation  (38). 
The  elastic  strain  increment  is  determined  by  equation  (31)  to  be 


dtE  .  o  .  IE-.?.1  -  t'ULLsJ. 

ai  u  Et+At  r 


The  total  strain  increment  is  thus 


dt  +  dt^ 


£  +  ftTKi  +  *aT 


It  can  be  seen  that  the  resulting  strain  increment  is  exactly  that 
predicted  by  equation  (37).  Therefore,  a  convergent  equilibrium 
Iteration  procedure  will  result  in  a  strain  and  stress  increment  in 
accordance  with  engineering  intuition.  Since  hardening  is  isotropic. 
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du  ■  0,  and  the  resulting  yield  surface  translation  tensor  is  v*. 

The  above  case  is  intended  to  illustrate  the  predictive  capability 
of  the  model  for  elastic  unloading. 

D.  Case  Number  ?  -  Temperature  Dependent  Mastic  Modulus 

Material  properties  and  load  history  for  the  second  sample  problem 
are  shown  in  Figure  7.  Hardening  is  assumed  to  be  kinematic.  The 
object  of  this  test  case  is  to  determine  whether  the  proper  stress 
Increment  is  predicted  when  the  elastic  modulus  is  temperature 
dependent.  As  an  initial  guess  the  strain  increment  is  estimated 
to  be. 

d.  »  1.2  «  ♦  a.\T  .  (41) 

y 

The  yield  surface  translation  at  the  start  of  the  step  is  found  to  be 

•  ,5  >>  .  The  loading  condition  [equation  (11]  will  be  satisfied 
y 

because  the  yield  surface  sire  is  decreasing  with  temperature.  Therefore, 
one  must  calculate  the  slope  of  the  uniavial  stress  plastic  strain  curve: 
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where  it  should  be  noted  that  3a/?T  is  calculated  using  the  equivalent 
plastic  strain  and  temperature  at  the  start  of  the  step.  Thus,  as  can 
be  seen  from  Figure  7,  the  model  predicts  the  correct  stress  increment 
when  the  elastic  modulus  is  temperature  dependent.  Since  hardening  is 
kinematic  equation  (34)  gives 


2(1.5  oy  -  .5  oy)  .  0  -  2  .  ov  .  gif  AT  , 
d“  =  0-5  -  .5  cy)  2  (1.5  oy  -  .5“T  *  3  .  (44) 

and  the  yield  surface  translation  increment  is  given  by  [equation  (35)]: 

1  °v 

da  =  3  (1.5  ay  -  .5  *y)  =  _  (45) 

A  graphical  representation  of  the  yield  surface  is  shown  in  Figure  C. 
C.  Case  Number  3  -  Verification  of  H' 

In  the  third  sample  case  the  bar  is  loaded  isothermally  to  some 
state  of  stress  outside  the  initial  yield  surface  and  then  slowly 
heated  and  mechanically  loaded,  as  shown  in  Figure  9.  Hardening  is 
assumed  to  be  isotropic.  The  elastic  modulus  is  independent  of  temperature 
so  that  this  case  checks  the  validity  of  H'.  From  Figure  9  one  can 
estimate  the  total  strain  increment  to  be 

de  =  4  Ey  +  aAT  .  (46) 

Since  equation  (32)  predicts  a  positive  stress  increment  and  the 
yield  surface  size  decreases  due  to  the  temperature  increase,  equation 
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(1)  predicts  loading.  Therefore,  it  Is  determined  that 


(-17) 


Since  the  elastic  modulus  does  not  change  the  second  and  third  terms 
disappear  in  equation  (33)  and  the  stress  increment  is  found  to  he 


da 


a 

.  3o  .  3o 

Ey  y  y 

i  i 
7  'y 

.  3e  .  3a  .  + 

>, 

>> 

-0 

,  .  -1 
y_.A 

•  3oy  7ZT 

AT  «  .25  o. 


3„y  *  .  3„y  .  3oy 


y 


(48) 


The  predicted  stress  increment  can  be  seen  to  be  in  accordance 
with  the  value  given  in  Figure  9.  Therefore,  H'  is  verified  for 
this  loading  case.  Since  the  hardening  is  isotropic  the  yield  surface 
translation  tensor  will  be  <|>,  as  can  be  seen  from  Figure  10. 

D.  Case  Number  4  -  Combined  Thermal  and  Mechanical  Loading 

The  final  sample  problem  illustrates  a  combined  loading  case  for 
an  axial  bar.  The  loading  conditions  and  material  properties  are 
shown  in  Figure  11.  Hardening  is  assumed  to  be  combined  isotropic- 
kinematic  with  n  *  .5.  From  the  figure  it  is  seen  that  the  strain 
Increment  can  be  estimated  to  be 


dc  ■  4.5  i  ♦  a AT.  (49) 

The  yield  surface  translation  at  the  start  of  the  step  is  determined  to 
be  «  ■  .25  o  .  An  evaluation  of  equation  (32)  will  once  again  predict 


loading  when  used  in  equation  (1).  Therefore,  H'  is  found  to  be 


The  stress  Increment,  from  equation  (33),  Is 


v  ^  .  2.5  o  .  2.5  o  .  y  ^ 

3c..  y  y  5 1.. 


2  i _ 3‘y  y 

V  y  y  ’  y 


J! _ 


2.5  ■  .  2.5  ,y 


■  2.5  oy  .  2.5  gy  ^  ^ 


a r?  •  2.5  o  .  2.5  «  +  v  ^  .  2.5  o  .  2.5  o  / 

,  y  y  y  J  cy  y  v 


1,5  ey) 


(3ty  '  1<5  ly} 


2.5  o 


y  *  I 


2.5  o. 


+  77 


2  °v  ?  °v 

5  .  2 . 5  o  .  2 . 5  o  +  ■*•  .  2 . 5  o  .  2 . 5  o 

9  cy  y  y  3  ey  y  y 


AT  =  .25  o. 


•  (51) 


Note  that  3o/DT  should  be  obtained  from  the  uniaxial  stress  plastic  strain 


diagram  rather  than  the  uniaxial  stress-strain  diagram.  This  can  be 


seen  by  noting  that  equation  (14)  Implies  that  do/dT  Is  evaluated  at 


constant  plastic  strain.  An  Inspection  of  Figure  11  shows  that  equation 


(51)  predicts  the  correct  stress  increment.  The  yield  surface  translation 


scalar  Is  [equation  (34)] 


du  -J2  (1.5  oy  -  .25  Oy)  .25  oy  -  2  .  1 .25  oy  .  ^  .AT 

-  2  .  1.25  ,y  .  ^  .  3.375  ryj^(1.5  „y  -  .25  oy)  2  .  .  (M) 


U 
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giving  a  yield  surface  translation  of 


d"  *  30  (1-5  °y  '  *25  °y)  *  24  °y  .  (53) 

A  careful  examination  of  the  given  data  will  show  that  the  yield  surface 
translation  tensor  guarantees  that  the  stress  tensor  remains  consistent 
with  the  yield  surface. 

CONCLUSION 

It  has  been  the  purpose  of  this  report  to  establish  the  procedure 
necessary  to  implement  the  temperature  dependent  elastic-plastic 
constitutive  law  previously  proposed  by  these  authors  into  a  finite 
element  code  for  structural  analysis.  In  so  doing,  we  have  presented 
some  modifications  to  the  previously  reported  law,  as  well  as  additional 
information  necessary  for  computer  applicability.  We  have  also  included 
a  computational  procedure,  and  we  have  shown  that  our  form  of  the 
constitutive  law  produces  correct,  results  for  several  uniaxial  sample 
cases.  In  every  sample  case  we  have  obtained  identical  results  in 
the  computer  code  AGGIE  I,  thus  verifying  not  only  the  constitutive  law 
but  the  finite  element  discretization  and  equilibrium  iteration  procedure 
as  well.  A  future  paper  will  compare  the  results  of  experiments  to 
analysis  for  specimens  subjected  to  cyclic  mechanical  and  thermal  loading. 
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FIGURE  3.  Construction  of  Yield  Surface  Size  Vs.  Uniaxial  Plastic 
Strain  Curves 


FIGURE  5.  Materials  Data  and  Load  History  for  Case  Number  1 
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FIGURE  9. 
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FIGURE  10.  Graphical  Representation  of  Yield  Surface  Expansion 
Case  Number  3 
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